unit relax;

interface
uses common;
Var
   anb	  : array[1..tdim,1..at_of_node] of boolean;
   n_dr	  : array[1..tdim,1..at_of_node] of v3; {node.atoms}
   bnn	  : array[1..nsites,1..nsites] of boolean;
   
procedure dparinit;
procedure relax_para_loc(nsteps: integer; xgamma,drmax:myreal;Var du : myreal);
procedure relax_para(nsteps : integer; xwt,xgamma,drmax:myreal;Var dr2,epot : myreal);
procedure cmc(nsteps,kk:integer; xwt,sigma:myreal; Var rate,epot,ueff,rmsd:myreal);
procedure relax(nsteps : integer; xgamma,drmax:myreal;Var epot : myreal);
procedure nntable(rplus	:myreal);
procedure uppw(xwt :myreal; Var epot: myreal; Var se:rsites; Var dse:rsites3);
procedure uppav(Var epot: myreal; Var se:rsites);
procedure uppaw(xwt :myreal; Var epot: myreal; Var se:rsites);
procedure elastic_E(xwt	: myreal;Var epot,elas:myreal;Var se,sel:rsites);

implementation
uses d04_header,ppot,etools,tools,outp;
procedure dparinit;
var i,j,o : integer;
begin
   for i:=1 to n_pxyz do begin
      for j:=1 to nn_pxyz[i] do begin ndps[i,j]:=0; d_pxyz[i,j]:=0.0; end;
   end;
end;
procedure nntable(rplus	:myreal);
Var
   i,j,o,m,n,ns	: integer;
   p1,p2,p3	: integer;
   r1,r2,v,w	: v3;
   p		: i3;
   bo		: boolean;
   ar2,rc2	: myreal;
   r		: rsites3;
begin	     
   {store pointers and atom-basis}
   ns:=0; for i:=1 to tnodes do for j:=1 to nan[i] do begin
      ns:=ns+1;
      para2xyz(n_para[i,j],0,1,0,r[ns]);
      for o:=1 to 3 do r[ns,o]:=bound(r[ns,o]+rt[i,o]);
   end;
   rc2:=(rcut+rplus)*(rcut+rplus);
   for i:=1 to ns do
      for j:=1 to ns do bnn[i,j]:=false;
   for i:=1 to ns do begin
      for j:=i to ns do begin
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  bo:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		  if not(bo) then begin
		     p[1]:=p1; p[2]:=p2; p[3]:=p3;
		     for o:=1 to 3 do v[o]:=bound(r[i,o]-r[j,o]);
		     for m:=1 to 3 do w[m]:=0.0;
		     for m:=1 to 3 do
			for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
		     ar2:=w[1]*w[1]+w[2]*w[2]+w[3]*w[3];
		     if (ar2<rc2) then begin
			bnn[i,j]:=true;
			bnn[j,i]:=true;
		     end;
		  end;
	       end;
      end;
   end;
end; { nntable }
procedure uppw(xwt :myreal; Var epot: myreal; Var se:rsites; Var dse:rsites3);
{returns energy, site energies, and partial derivatives in cartezian basis}
Var
   i,j,k,m,n,o	: integer;
   n1,n2,a1,a2	: integer;
   nat,p1,p2,p3	: integer;
   io1,io2	: integer;
   itp1,itp2	: integer;
   w,v		: v3;
   p		: i3;
   a		: array[1..nsites] of i2;
   b		: boolean;
   q		: array[1..nsites] of integer;
   r		: rsites3;
   ar,u,ud,udd	: myreal;
   u0,ud0,udd0	: myreal;
   wt,xw	: myreal;
begin
   nat:=0; k:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do begin
	 k:=k+1;
	 if (n_a[i,j,1]>0) then begin
	    nat:=nat+1; a[nat,1]:=i; a[nat,2]:=j;
	    para2xyz(n_para[i,j],n_a[i,j,1],1,1,n_r[i,j]);
	    for o:=1 to 3 do r[nat,o]:=bound(n_r[i,j,o]+rt[i,o]);
	    q[nat]:=k;
	 end;
      end;
   epot:=0.0;
   for i:=1 to nat do begin
      se[i]:=0.0;
      for o:=1 to 3 do dse[i,o]:=0.0;
   end;
   for i:=1 to nat do begin
      n1:=a[i,1]; a1:=a[i,2]; itp1:=n_a[n1,a1,1]; io1:=n_a[n1,a1,2];
      for j:=i to nat do if bnn[q[i],q[j]] then begin
	 n2:=a[j,1]; a2:=a[j,2]; itp2:=n_a[n2,a2,1]; io2:=n_a[n2,a2,2];
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  b:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		  if not(b) then begin
		     for o:=1 to 3 do v[o]:=bound(r[j,o]-r[i,o]);
		     p[1]:=p1; p[2]:=p2; p[3]:=p3;
		     for m:=1 to 3 do w[m]:=0.0;
		     for m:=1 to 3 do
			for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
		     ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
		     if (ar<rcut) then begin
			upp(itp1,itp2,ar,u0,ud0,udd0);
			wt:=2/(naoch[io1,itp1]+naoch[io2,itp2]);
			xw:=(1+(wt-1)*xwt);
			u:=u0*xw;
			ud:=ud0*xw;
			udd:=udd0*xw;
			for o:=1 to 3 do begin  {cartezian partial derivatives}
			   dse[i,o]:=dse[i,o]+w[o]*ud;
			   dse[j,o]:=dse[j,o]-w[o]*ud;
			end;
			epot:=epot+u;
			if (i=j) then epot:=epot-u/2.0;
			se[i]:=se[i]+u/2.0;
			if (i<>j) then se[j]:=se[j]+u/2.0;
		     end;
		  end;
	       end;
      end;
   end;
end; { uppa }
procedure elastic_E(xwt	: myreal;Var epot,elas:myreal;Var se,sel:rsites);
{energy datapoint; sets positional derivatives from the pair potentials}
{nsteps xwt xgamma drmax dr2 epot}
{ASSUMES BNN table is initialized!}
Var
   i,j,k,l,o,m,n	   : integer;
   bbij,cmp,c,els	   : myreal;
   io,jo,nat,io1,io2	   : integer;
   p1,p2,p3		   : integer;
   p			   : i3;
   q			   : array[1..nsites] of integer;
   ar,ar2,u,ud,udd	   : myreal;
   u0,ud0,udd0		   : myreal;
   v,w,w1,drt,dpr,de,ri,rj : v3;
   bo			   : boolean;
   a			   : array[1..nsites] of i2;
   bb			   : array[1..nsites] of v33;
   br			   : v33;
   b			   : rsites3;
   itp1,itp2,n1,n2,a1,a2   : integer;
   wt,xw		   : myreal;
begin
   {this is how the mass come in: bbij*srm(ita)*srm(jta)}
   {store pointers and atom-basis}
   nat:=0; k:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do begin
	 k:=k+1;
	 if (n_a[i,j,1]>0) then begin
	    nat:=nat+1; a[nat,1]:=i; a[nat,2]:=j;
	    para2xyz(n_para[i,j],n_a[i,j,1],1,1,n_r[i,j]);
	    q[nat]:=k;
	 end;
      end;
   epot:=0.0;
   for i:=1 to nat do begin
      se[i]:=0;
      for io:=1 to 3 do begin
	 b[i,io]:=0.0;
	 for jo:=1 to 3 do bb[i,io,jo]:=0.0;
      end;
   end;
   dparinit;
   for i:=1 to nat do begin
      n1:=a[i,1]; a1:=a[i,2]; itp1:=n_a[n1,a1,1]; io1:=n_a[n1,a1,2];
      for j:=i to nat do if bnn[q[i],q[j]] then begin 
	 n2:=a[j,1]; a2:=a[j,2]; itp2:=n_a[n2,a2,1]; io2:=n_a[n2,a2,2];
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  bo:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		  if not(bo) then begin
		     p[1]:=p1; p[2]:=p2; p[3]:=p3;
		     for o:=1 to 3 do v[o]:=bound(n_r[n1,a1,o]-n_r[n2,a2,o]+rt[n1,o]-rt[n2,o]);
		     for m:=1 to 3 do w[m]:=0.0;
		     for m:=1 to 3 do
			for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
		     ar2:=w[1]*w[1]+w[2]*w[2]+w[3]*w[3];
		     ar:=sqrt(ar2);
		     if (ar<rcut) and (ar>0.001) then begin
			{wt:=2/(naoch[io1,itp1]+naoch[io2,itp2]);}
			upp(itp1,itp2,ar,u0,ud0,udd0);
			wt:=2/(naoch[io1,itp1]+naoch[io2,itp2]);
			xw:=(1+(wt-1)*xwt);
			u:=u0*xw;
			ud:=ud0*xw;
			udd:=udd0*xw;
			se[i]:=se[i]+u/2.0;
			if (i<>j) then se[j]:=se[j]+u/2.0;
			for io:=1 to 3 do begin
			   b[i,io]:=b[i,io]+ud*w[io];
			   for jo:=1 to 3 do begin
			      bbij:=(ud-udd)*w[io]*w[jo]/ar2;
			      IF (IO=JO) then BBIJ:=BBIJ-ud;
			      BB[i,IO,JO]:=BB[i,IO,JO]+BBIJ;
			   end;
			end;
			for o:=1 to 3 do w[o]:=-w[o];
			for io:=1 to 3 do begin
			   b[j,io]:=b[j,io]+ud*w[io];
			   for jo:=1 to 3 do begin
			      bbij:=(ud-udd)*w[io]*w[jo]/ar2;
			      IF (IO=JO) then BBIJ:=BBIJ-ud;
			      BB[j,IO,JO]:=BB[j,IO,JO]+BBIJ;
			   end;
			end;
		     end;
		  end;
	       end;
      end;
   end;
   elas:=0.0;
   for i:=1 to nat do begin
      p1:=a[i,1]; p2:=a[i,2];
      {elstic displacement; bb_ij is K_ij and b is force}
      mxle3(bb[i],b[i],drt);
      sel[i]:=0.0;
      for io:=1 to 3 do
	 for jo:=1 to 3 do
	    sel[i]:=sel[i]+drt[io]*bb[i,io,jo]*drt[jo];
      sel[i]:=-sel[i]/2;
      elas:=elas+sel[i];
      epot:=epot+se[i];
      c:=sqrt(abs(sel[i])); {square root of the mean force constant}
      for io:=1 to 3 do drt[io]:=drt[io]*c;
      get_atom_basis(n_para[p1,p2],br);
      mxle3(br,drt,dpr); {parameter space basis}
      add2para(n_para[p1,p2],n_a[p1,p2,1],dpr);
   end;
end; { elastic_E }
procedure uppav(Var epot: myreal; Var se:rsites);
{uppa-v, "v" for energy only, not 1st and 2nd derivatives, to save time}
{returns energy, and site energies, no derivatives}
Var
   i,j,k,m,n,o	: integer;
   n1,n2,a1,a2	: integer;
   nat,p1,p2,p3	: integer;
   itp1,itp2	: integer;
   w,v		: v3;
   p		: i3;
   a		: array[1..nsites] of i2;
   b		: boolean;
   q		: array[1..nsites] of integer;
   r		: rsites3;
   ar,u,ud,udd	: myreal;
begin
   nat:=0; k:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do begin
	 k:=k+1;
	 if (n_a[i,j,1]>0) then begin
	    nat:=nat+1; a[nat,1]:=i; a[nat,2]:=j;
	    para2xyz(n_para[i,j],n_a[i,j,1],1,1,n_r[i,j]);
	    for o:=1 to 3 do r[nat,o]:=bound(n_r[i,j,o]+rt[i,o]);
	    q[nat]:=k;
	 end;
      end;
   epot:=0.0;
   for i:=1 to nat do se[i]:=0.0;
   for i:=1 to nat do begin
      n1:=a[i,1]; a1:=a[i,2]; itp1:=n_a[n1,a1,1];
      for j:=i to nat do if bnn[q[i],q[j]] then begin
	 n2:=a[j,1]; a2:=a[j,2]; itp2:=n_a[n2,a2,1];
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  b:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		  if not(b) then begin
		     for o:=1 to 3 do v[o]:=bound(r[j,o]-r[i,o]);
		     p[1]:=p1; p[2]:=p2; p[3]:=p3;
		     for m:=1 to 3 do w[m]:=0.0;
		     for m:=1 to 3 do
			for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
		     ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
		     if (ar<rcut) then begin
			uppv(itp1,itp2,ar,u);
			epot:=epot+u;
			if (i=j) then epot:=epot-u/2.0;
			se[i]:=se[i]+u/2.0;
			if (i<>j) then se[j]:=se[j]+u/2.0;
		     end;
		  end;
	       end;
      end;
   end;
end; { uppav }
procedure uppaw(xwt :myreal; Var epot: myreal; Var se:rsites);
{uppa-w, "w" for WEIGHTED, using weights naoch for each pair contributing to the energy}
{returns energy, and site energies, no derivatives}
Var
   i,j,k,m,n,o	: integer;
   n1,n2,a1,a2	: integer;
   nat,p1,p2,p3	: integer;
   itp1,itp2	: integer;
   io1,io2	: integer;
   w,v		: v3;
   p		: i3;
   a		: array[1..nsites] of i2;
   b		: boolean;
   q		: array[1..nsites] of integer;
   r		: rsites3;
   ar,u,u0,wt	: myreal;
begin
   nat:=0; k:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do begin
	 k:=k+1;
	 if (n_a[i,j,1]>0) then begin
	    nat:=nat+1; a[nat,1]:=i; a[nat,2]:=j;
	    para2xyz(n_para[i,j],n_a[i,j,1],1,1,n_r[i,j]);
	    for o:=1 to 3 do r[nat,o]:=bound(n_r[i,j,o]+rt[i,o]);
	    q[nat]:=k;
	 end;
      end;
   epot:=0.0;
   for i:=1 to nat do se[i]:=0.0;
   for i:=1 to nat do begin
      n1:=a[i,1]; a1:=a[i,2]; itp1:=n_a[n1,a1,1]; io1:=n_a[n1,a1,2];
      for j:=i to nat do if bnn[q[i],q[j]] then begin
	 n2:=a[j,1]; a2:=a[j,2]; itp2:=n_a[n2,a2,1]; io2:=n_a[n2,a2,2];
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  b:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		  if not(b) then begin
		     for o:=1 to 3 do v[o]:=bound(r[j,o]-r[i,o]);
		     p[1]:=p1; p[2]:=p2; p[3]:=p3;
		     for m:=1 to 3 do w[m]:=0.0;
		     for m:=1 to 3 do
			for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
		     ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
		     if (ar<rcut) then begin
			wt:=2/(naoch[io1,itp1]+naoch[io2,itp2]);
			uppv(itp1,itp2,ar,u0);
			u:=u0*(1+(wt-1)*xwt);
			epot:=epot+u;
			if (i=j) then epot:=epot-u/2.0;
			se[i]:=se[i]+u/2.0;
			if (i<>j) then se[j]:=se[j]+u/2.0;
		     end;
		  end;
	       end;
      end;
   end;
end; { uppav }
procedure relax_para_loc(nsteps:integer; xgamma,drmax:myreal;Var du : myreal);
Var
   i,j,k,l,o,m,n,il	     : integer;
   bbij,ep0,epot	     : myreal;
   io,jo,nat		     : integer;
   i1,i2,p1,p2,p3	     : integer;
   p			     : i3;
   ar,ar2,u,ud,udd	     : myreal;
   cf,drm2		     : myreal;
   v,w,w1,drt,drx,ri,rj	     : v3;
   bo			     : boolean;
   a			     : array[1..nsites,1..3] of integer;
   bb			     : array[1..nsites] of v33;
   br			     : v33;
   b,dpr,dra		     : rsites3;
   itp1,itp2,n1,n2,a1,a2,nre : integer;
   oa			     : array[1..t_orb] of boolean;
begin
   for i:=1 to t_orb do oa[i]:=false;
   drm2:=drmax*drmax;
   {store pointers and atom-basis}
   nat:=0; nre:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
	 nat:=nat+1; a[nat,1]:=i; a[nat,2]:=j;
	 if (o_rel[n_a[i,j,2]]) then a[nat,3]:=0 else begin
	    nre:=nre+1;
	    a[nat,3]:=1; {relax status}
	    oa[n_a[i,j,2]]:=true;
	 end;
      end;
   {writeln('relax: ',nat:1,' atoms.');}
   if (nre>0) then begin
      if debug[8] then begin
	 write('LOC_REL: adding orbits ');
	 for i:=1 to t_orb do if oa[i] then write(i:1,' ');
	 writeln(' -> ',nre:1,' atoms');
      end;
      for il:=1 to nsteps do begin
	 epot:=0.0;
	 for i:=1 to nat do begin
	    for io:=1 to 3 do begin
	       b[i,io]:=0.0;
	       for jo:=1 to 3 do bb[i,io,jo]:=0.0;
	    end;
	 end;
	 dparinit;
	 for i:=1 to nat do if (a[i,3]=1) then begin
	    for j:=i to nat do begin
	       for p1:=-iext[1] to iext[1] do
		  for p2:=-iext[2] to iext[2] do
		     for p3:=-iext[3] to iext[3] do begin
			n1:=a[i,1]; a1:=a[i,2];
			n2:=a[j,1]; a2:=a[j,2];
			itp1:=n_a[n1,a1,1];
			itp2:=n_a[n2,a2,1];
			bo:=((p1=0) and (p2=0) and (p3=0) and (i=j));
			if not(bo) then begin
			   p[1]:=p1; p[2]:=p2; p[3]:=p3;
			   for o:=1 to 3 do v[o]:=bound(n_r[n1,a1,o]-n_r[n2,a2,o]+rt[n1,o]-rt[n2,o]);
			   for m:=1 to 3 do w[m]:=0.0;
			   for m:=1 to 3 do
			      for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
			   ar2:=w[1]*w[1]+w[2]*w[2]+w[3]*w[3];
			   ar:=sqrt(ar2);
			   if (ar<rcut) and (ar>0.01) then begin
			      upp(itp1,itp2,ar,u,ud,udd);
			      epot:=epot+u;
			      if (i=j) then epot:=epot-u/2.0;
			      for io:=1 to 3 do begin
				 {force: notice ud is already divided by R}
				 b[i,io]:=b[i,io]+ud*w[io]; 
				 for jo:=1 to 3 do begin
				    bbij:=(ud-udd)*w[io]*w[jo]/ar2;
				    IF (IO=JO) then BBIJ:=BBIJ-ud;
				    BB[i,IO,JO]:=BB[i,IO,JO]+BBIJ;
				 end;
			      end;
			      if (a[j,3]=1) then begin
				 for o:=1 to 3 do w[o]:=-w[o];
				 for io:=1 to 3 do begin
				    b[j,io]:=b[j,io]+ud*w[io];
				    for jo:=1 to 3 do begin
				       bbij:=(ud-udd)*w[io]*w[jo]/ar2;
				       IF (IO=JO) then BBIJ:=BBIJ-ud;
				       BB[j,IO,JO]:=BB[j,IO,JO]+BBIJ;
				    end;
				 end;
			      end;
			   end;
			end;
		     end;
	    end;
	 end;
	 if (il=1) then ep0:=epot;
	 {move}
	 if (il<nsteps) then begin {last step only calculates energy}
	    for i:=1 to nat do if (a[i,3]=1) then begin
	       p1:=a[i,1]; p2:=a[i,2];
	       mxle3(bb[i],b[i],drt);
	       {check if displacement isnt very large; drt is in A}
	       ar2:=0.0; for o:=1 to 3 do ar2:=ar2+drt[o]*drt[o];
	       cf:=1.0;
	       if (ar2>drm2) then begin
		  cf:=sqrt(drm2/ar2);
		  for o:=1 to 3 do drt[o]:=drt[o]*cf;
	       end;
	       for o:=1 to 3 do drt[o]:=drt[o]*xgamma;
	       {for o:=1 to 6 do write(n_para[p1,p2,o]:1,' '); writeln;}
	       get_atom_basis(n_para[p1,p2],br);
	       mxle3(br,drt,dpr[i]); {parameters}
	       mxle3(bl3,drt,dra[i]); {crystallographic}
	       {write(i:1,' ',n_a[p1,p2,1]:1,' ');
	       for o:=1 to 3 do write(dpr[i,o]:1:5,' '); writeln;}
	       if (o_rel[n_a[p1,p2,2]]) then add2para(n_para[p1,p2],n_a[p1,p2,1],dpr[i]);
	    end;
	    {update parameters}
	    for i:=1 to n_pxyz do
	       for j:=1 to nn_pxyz[i] do if (ndps[i,j]>0) then begin
		  d_pxyz[i,j]:=d_pxyz[i,j]/ndps[i,j];
		  {writeln(i:1,' ',j:1,' ',ndps[i,j]:1,' ',dpxyz[i,j]:1:5,' ',d_pxyz[i,j]:1:5);}
		  dpxyz[i,j]:=dpxyz[i,j]+d_pxyz[i,j];
	       end;
	    ipos:=1; iposd:=1;
	    for i:=1 to nat do if (a[i,3]=1) then begin {update coordinates}
	       p1:=a[i,1]; p2:=a[i,2];
	       o_rel[n_a[p1,p2,2]]:=true;
	       para2xyz(n_para[p1,p2],n_a[p1,p2,1],ipos,iposd,n_r[p1,p2]);
	    end;
	 end;
	 if debug[8] then begin
	    writeln('LOC_REL ',il:2,'  dU ',(epot-ep0):1:6);
	 end;
      end; {main nsteps loop}
      {UPPN(epot);} {update energy after the last step}
      du:=du+epot-ep0;
   end;
end; { relax_para_loc }
procedure cmc(nsteps,kk:integer; xwt,sigma:myreal; Var rate,epot,ueff,rmsd:myreal);
{ASSUMES BNN table is initialized!}
{displace only one orbit at a time}
Var
   i,j,k,l,o,il,norb : integer;
   ep0,ep1,eini	     : myreal;
   p1,p2	     : integer;
   cf,cmp,sg2,cs2,cs : myreal;
   du0,x,y	     : myreal;
   v,w,w1,drt,dpr    : v3;
   br		     : v33;
   se		     : rsites;
   accept	     : boolean;
   nacc,lab	     : integer;
   m1,m2,nao,naao    : integer;
procedure update_par(m1,m2,isg: integer;dr:v3);
Var
   o,k,pp : integer;
begin	  
   for o:=1 to 3 do if (n_para[m1,m2,o]>0) then begin
      pp:=n_para[m1,m2,o];
      for k:=1 to nn_pxyz[pp] do if (n_a[m1,m2,1]=ch_pxyz[pp,k]) then begin
	 dpxyz[pp,k]:=dpxyz[pp,k]+isg*dr[o];
      end;
   end;
end;
begin
   lab:=9;
   eini:=epot;
   {writeln('eini ',eini:1:6);}
   och_status;
   nacc:=0;rate:=0.0;naao:=0;
   sg2:=2*sigma*sigma;
   if debug[8] then writeln(sname[9],': mcs beta ',nsteps:1,' ',beta[kk]:1:4);
   cs:=0.0; cs2:=0.0;
   UPPAW(xwt,epot,se);
   for il:=1 to nsteps*natoms do begin
      ep0:=epot;
      repeat
	 cmp:=randu(seed)*tnodes; {pick randomly tiling node}
	 m1:=round(cmp+0.5); 
	 cmp:=randu(seed)*nan[m1]; {pick randomly an atom}
	 m2:=round(cmp+0.5);
      until n_a[m1,m2,1]>0;
      nao:=naoch[n_a[m1,m2,2],n_a[m1,m2,1]];
      for o:=1 to 3 do drt[o]:=randu(seed)-0.5;
      repeat {gaussian random length (0,drmax)}
	 cf:=randu(seed)*sigma*5; {should be big enough}
	 y:=randu(seed);
      until (exp(-cf*cf/sigma/sigma/2)>y);
      {cf:=cf/naoch[n_a[m1,m2,2],n_a[m1,m2,1]];}
      cmp:=absv3(drt);
      {writeln('cf ',cf:1:4,' ',cmp:1:4);}
      for o:=1 to 3 do drt[o]:=drt[o]*cf/cmp;
      get_atom_basis(n_para[m1,m2],br);
      mxle3(br,drt,dpr); {parameters basis}
      update_par(m1,m2,1,dpr);
      {energy}
      UPPAW(xwt,ep1,se);
      du0:=(ep1-ep0)*beta[kk];
      accept:=true;
      if(du0>0)then begin
	 x:=exp(-du0);
	 y:=randu(seed);
	 if (y>x) then begin {reject}
	    update_par(m1,m2,-1,dpr);
	    accept:=false;
	 end else begin {accept}
	 end;
      end;
      if(accept)then begin
	 if ((il mod save_e)=0) then begin
	    write(recu,(ep1-ep0):1:6,' ');
	    writeln(recu,'   ',kk:1,' ',lab:1);
	 end;
	 nacc:=nacc+1;
	 epot:=ep1;
	 cs2:=cs2+cf*cf*nao;
	 naao:=naao+nao;
	 mystat(kk,lab,epot);
	 {if(epot<emin)then begin
	    export_config_if(emin,epot/natoms,beta[kk],sname[9]);
	 end;}
      end;
      if (debug[8] and accept) then writeln(il:2,'   ',epot/natoms:1:6,' ',(ep1-ep0):1:4);
   end; {main nsteps loop}
   ueff:=epot/natoms;
   rmsd:=sqrt(cs2/naao);
   UPPAV(epot,se); {update energy after the last step}
   if debug[8] then begin
      writeln('dU(tot) ',(epot-eini):1:6,'  dU/at ',(epot-eini)/natoms:1:6);
      if(naao>0)then writeln('rms-displacement: ',sqrt(cs2/naao):1:4);
   end;
   rate:=nacc/nsteps/natoms;
end; { cmc }
procedure relax_para(nsteps : integer; xwt,xgamma,drmax:myreal;Var dr2,epot : myreal);
{ASSUMES BNN table is initialized!}
{One attempt is for all orbits}
Var
   i,j,i0,j0,k,l,o,m,n,il : integer;
   bbij,ep0,ep1,cmp	  : myreal;
   io,jo,nat,io1,io2	  : integer;
   i1,i2,p1,p2,p3,nb	  : integer;
   p			  : i3;
   q			  : array[1..nsites] of integer;
   ar,ar2,u,ud,udd	  : myreal;
   u0,ud0,udd0		  : myreal;
   cf,drm2		  : myreal;
   v,w,w1,drt,drx,ri,rj	  : v3;
   bo			  : boolean;
   a			  : array[1..nsites,1..2] of integer;
   bb			  : array[1..nsites] of v33;
   br			  : v33;
   b,dpr,dra		  : rsites3;
   itp1,itp2,n1,n2,a1,a2  : integer;
   wt,xw		  : myreal;
   se			  : rsites;
begin
   drm2:=drmax*drmax;
   dr2:=0.0;
   {store pointers and atom-basis}
   nat:=0; k:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do begin
	 k:=k+1;
	 if (n_a[i,j,1]>0) then begin
	    nat:=nat+1; a[nat,1]:=i; a[nat,2]:=j;
	    q[nat]:=k;
	 end;
      end;
   {writeln('relax: ',nat:1,' atoms, rcplus=',rcplus:1:1);}
   for il:=1 to nsteps+1 do begin
      ep1:=0.0;
      for i:=1 to nat do begin
	 for io:=1 to 3 do begin
	    b[i,io]:=0.0;
	    for jo:=1 to 3 do bb[i,io,jo]:=0.0;
	 end;
      end;
      dparinit;
      nb:=0;
      for i:=1 to nat do begin
	 n1:=a[i,1]; a1:=a[i,2]; itp1:=n_a[n1,a1,1]; io1:=n_a[n1,a1,2];
	 for j:=i to nat do if bnn[q[i],q[j]] then begin 
	    n2:=a[j,1]; a2:=a[j,2]; itp2:=n_a[n2,a2,1]; io2:=n_a[n2,a2,2];
	    for p1:=-iext[1] to iext[1] do
 for p2:=-iext[2] to iext[2] do
		  for p3:=-iext[3] to iext[3] do begin
		     bo:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		     if not(bo) then begin
			p[1]:=p1; p[2]:=p2; p[3]:=p3;
			for o:=1 to 3 do v[o]:=bound(n_r[n1,a1,o]-n_r[n2,a2,o]+rt[n1,o]-rt[n2,o]);
			for m:=1 to 3 do w[m]:=0.0;
			for m:=1 to 3 do
			   for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
			ar2:=w[1]*w[1]+w[2]*w[2]+w[3]*w[3];
			ar:=sqrt(ar2);
			if (ar<rcut) and (ar>0.001) then begin
			   nb:=nb+1;
			   {wt:=2/(naoch[io1,itp1]+naoch[io2,itp2]);}
			   upp(itp1,itp2,ar,u0,ud0,udd0);
			   wt:=2/(naoch[io1,itp1]+naoch[io2,itp2]);
			   xw:=(1+(wt-1)*xwt);
			   u:=u0*xw;
			   ud:=ud0*xw;
			   udd:=udd0*xw;
			   ep1:=ep1+u;
			   if(i=j) then ep1:=ep1-u/2;
			   for io:=1 to 3 do begin
			      b[i,io]:=b[i,io]+ud*w[io];
			      for jo:=1 to 3 do begin
				 bbij:=(ud-udd)*w[io]*w[jo]/ar2;
				 IF (IO=JO) then BBIJ:=BBIJ-ud;
				 BB[i,IO,JO]:=BB[i,IO,JO]+BBIJ;
			      end;
			   end;
			   for o:=1 to 3 do w[o]:=-w[o];
			   for io:=1 to 3 do begin
			      b[j,io]:=b[j,io]+ud*w[io];
			      for jo:=1 to 3 do begin
				 bbij:=(ud-udd)*w[io]*w[jo]/ar2;
				 IF (IO=JO) then BBIJ:=BBIJ-ud;
				 BB[j,IO,JO]:=BB[j,IO,JO]+BBIJ;
			      end;
			   end;
			end;
		     end;
		  end;
	 end;
      end;
      if (il=1) then ep0:=ep1;
      {move}
      if(il<=nsteps) then begin
	 for i:=1 to nat do begin
	    p1:=a[i,1]; p2:=a[i,2];
	    {write(p1:1,' ',p2:1,'  '); for o:=1 to 3 do write(b[i,o]:1:3,' '); writeln;}
	    mxle3(bb[i],b[i],drt);
	    for o:=1 to 3 do drt[o]:=drt[o]*xgamma;
	    {check if displacement isnt very large; drt is in A}
	    ar2:=0.0; for o:=1 to 3 do ar2:=ar2+drt[o]*drt[o];
	    cf:=1.0;
	    if (ar2>drm2) then begin
	       cf:=sqrt(drm2/ar2);
	       for o:=1 to 3 do drt[o]:=drt[o]*cf;
	    end;
	    cmp:=absv3(drt); dr2:=dr2+cmp*cmp;
	    {write(p1:1,' ',p2:1,'  '); for o:=1 to 6 do write(n_para[p1,p2,o]:1,' '); writeln;}
	    get_atom_basis(n_para[p1,p2],br);
	    mxle3(br,drt,dpr[i]); {parameters}
	    {mxle3(bl3,drt,dra[i]);} {crystallographic}
	    {write(i:1,' ',n_a[p1,p2,1]:1,' ');
	     for o:=1 to 3 do write(dpr[i,o]:1:5,' '); writeln;}
	    if (o_rel[n_a[p1,p2,2]]) then add2para(n_para[p1,p2],n_a[p1,p2,1],dpr[i]);
	 end;
	 if not(swap_man) then begin
	    {update parameters}
	    for i:=1 to n_pxyz do
	       for j:=1 to nn_pxyz[i] do if (ndps[i,j]>0) then begin
		  d_pxyz[i,j]:=d_pxyz[i,j]/ndps[i,j];
		  {writeln(i:1,' ',j:1,' ',ndps[i,j]:1,' ',dpxyz[i,j]:1:5,' ',d_pxyz[i,j]:1:5);}
		  dpxyz[i,j]:=dpxyz[i,j]+d_pxyz[i,j];
	       end;
	 end;
	 for i:=1 to nat do begin {update coordinates}
	    p1:=a[i,1]; p2:=a[i,2];
	    para2xyz(n_para[p1,p2],n_a[p1,p2,1],1,1,n_r[p1,p2]);
	 end;
	 if debug[8] then writeln(il:2,'   ',ep1/nat:1:6);
      end; {main nsteps loop}
   end;
   if debug[8] then writeln('dU(tot) ',(ep1-ep0):1:6,'  dU/at ',(ep1-ep0)/nat:1:6);
   xw:=(1+(nat/maxorb-1)*xwt);
   epot:=ep1;
   dr2:=dr2/nsteps;
end; { relax_para }
procedure relax(nsteps : integer; xgamma,drmax:myreal;Var epot : myreal);
Var
   i,j,k,l,o,m,n,il	 : integer;
   bbij,ep0		 : myreal;
   io,jo,nat		 : integer;
   i1,i2,p1,p2,p3	 : integer;
   p			 : i3;
   ar,ar2,u,ud,udd	 : myreal;
   cf,drm2		 : myreal;
   v,w,w1,drt,drx,ri,rj	 : v3;
   bo			 : boolean;
   a			 : array[1..nsites,1..2] of integer;
   bb			 : array[1..nsites] of v33;
   br			 : v33;
   b,dpr,dra		 : rsites3;
   se			 : rsites;
   itp1,itp2,n1,n2,a1,a2 : integer;
begin
   ep0:=epot;
   drm2:=drmax*drmax;
   {store pointers}
   nat:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
	 nat:=nat+1; a[nat,1]:=i; a[nat,2]:=j;
      end;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do
	 for o:=1 to 3 do n_dr[i,j,o]:=0.0;
   {writeln('relax: ',nat:1,' atoms.');}
   for il:=1 to nsteps do begin
      epot:=0.0;
      for i:=1 to nat do begin
	 for io:=1 to 3 do begin
	    b[i,io]:=0.0;
	    for jo:=1 to 3 do bb[i,io,jo]:=0.0;
	 end;
      end;
      for i:=1 to nat do begin
	 for j:=i to nat do begin
	    for p1:=-iext[1] to iext[1] do
	       for p2:=-iext[2] to iext[2] do
		  for p3:=-iext[3] to iext[3] do begin
		     n1:=a[i,1]; a1:=a[i,2];
		     n2:=a[j,1]; a2:=a[j,2];
		     itp1:=n_a[n1,a1,1];
		     itp2:=n_a[n2,a2,1];
		     bo:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		     if not(bo) then begin
			p[1]:=p1; p[2]:=p2; p[3]:=p3;
			for o:=1 to 3 do v[o]:=bound(n_r[n1,a1,o]-n_r[n2,a2,o]+rt[n1,o]-rt[n2,o]);
			for m:=1 to 3 do w[m]:=0.0;
			for m:=1 to 3 do
			   for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl3[m,n];
			ar2:=w[1]*w[1]+w[2]*w[2]+w[3]*w[3];
			ar:=sqrt(ar2);
			if (ar<rcut) and (ar>0.01) then begin
			   upp(itp1,itp2,ar,u,ud,udd);
			   epot:=epot+u;
			   if (i=j) then epot:=epot-u/2.0;
			   for io:=1 to 3 do begin
			      b[i,io]:=b[i,io]+ud*w[io];
			      for jo:=1 to 3 do begin
				 bbij:=(ud-udd)*w[io]*w[jo]/ar2;
				 IF (IO=JO) then BBIJ:=BBIJ-ud;
				 BB[i,IO,JO]:=BB[i,IO,JO]+BBIJ;
			      end;
			   end;
			   for o:=1 to 3 do w[o]:=-w[o];
			   for io:=1 to 3 do begin
			      b[j,io]:=b[j,io]+ud*w[io];
			      for jo:=1 to 3 do begin
				 bbij:=(ud-udd)*w[io]*w[jo]/ar2;
				 IF (IO=JO) then BBIJ:=BBIJ-ud;
				 BB[j,IO,JO]:=BB[j,IO,JO]+BBIJ;
			      end;
			   end;
			end;
		     end;
		  end;
	 end;
      end;
      for i:=1 to nat do begin {MOVE atoms}
	 p1:=a[i,1]; p2:=a[i,2];
	 mxle3(bb[i],b[i],drt);
	 {check if displacement isnt very large; drt is in A}
	 ar2:=0.0; for o:=1 to 3 do ar2:=ar2+drt[o]*drt[o];
	 cf:=1.0;
	 if (ar2>drm2) then begin
	    cf:=sqrt(drm2/ar2);
	    for o:=1 to 3 do drt[o]:=drt[o]*cf;
	 end;
	 for o:=1 to 3 do drt[o]:=drt[o]*xgamma;
	 mxle3(bl3,drt,dra[i]); {crystallographic}
	 if anb[p1,p2] then for o:=1 to 3 do begin
	    n_r[p1,p2,o]:=n_r[p1,p2,o]+dra[i,o];
	    n_dr[p1,p2,o]:=n_dr[p1,p2,o]+dra[i,o];
	 end;
      end;
      if debug[8] then writeln(il:2,'  ',epot/nat:1:6);
   end; {main nsteps loop}
   UPPAV(epot,se); {update energy after the last step}
end; { relax }

end.
